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Abstract 

We generalize a recently developed method for accelerated Monte Carlo calculation 
of path integrals to the physically relevant case of generic many-body systems. This 
is done by developing an analytic procedure for constructing a hierarchy of effective 
actions leading to improvements in convergence of A^-fold discretized many-body 
path integral expressions from to l/N^ for generic p. In this paper we present 
explicit solutions within this hierarchy up to level p = 5. Using this we calculate the 
low lying energy levels of a two particle model with quartic interactions for several 
values of coupling and demonstrate agreement with analytical results governing 
the increase in efficiency of the new method. The applicability of the developed 
scheme is further extended to the calculation of energy expectation values through 
the construction of associated energy estimators exhibiting the same speedup in 
convergence. 

Key words: Path integral, Effective action, Many-body system, Energy estimator 
PACS: 05.30.-d, 03.65.Db, 05.10.Ln 



1 Introduction 



Originally introduced in quantum mechanics p]l2] and later most widely used 
in high energy theory [3M] and condensed matter physics [51I6] . path integrals 
have become important tools throughout the physical sciences from atomic, 
molecular and nuclear physics, to chemistry and biophysics. Moreover, path 
integrals are starting to play important roles in several areas of mathematics 
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and in modern finance [7], especially in option pricing applications p[9] . The 
downside of the formalism is that the mathematical properties of path integrals 
are not sufficiently well understood, and that an extremely small number of 
path integrals can be solved exactly [TD]. An extensive overview of the path 
integral formalism and its various applications can be found in 

Beside their key position in analytical approaches to quantum theory, path 
integrals have an important role in direct numerical simulations of realistic 
many-body systems. Such simulations have made possible practical compar- 
isons of the predictions of various theoretical models with the results of asso- 
ciated experiments. Further, they have led to a deeper understanding of the 
complex physical phenomena involved [12] . 

The starting point in all such calculations is the time-sliced expression for the 
general quantum-mechanical transition amplitude p], 

AN{a,b;T) = ^-jjj^ f dqi ■ ■ ■ dq^ -i e~^^. (1) 



from initial state \a) to final state for time interval T, where N is the 
number of time slices, e = T/N and is the naively discretized action 
for a system of M non-relativistic particles in d spatial dimensions. The 
N oo limit of the above discretized amplitude gives the continuum am- 
plitude A{a,b;T) = (6|e~'^"^|a). The evaluation of these types of discretized 
expressions is handled by a variety of well developed numerical integration 
methods, and this is one of the principle reasons for the popularity of the 
path integral approach in numerical simulations. 

The performance of numerical algorithms for the calculation of path integrals 
is determined by the efficiency of the applied integration techniques, as well as 
by the speed of convergence of the discretized amplitudes to the continuum. 
The first aspect has been successfully dealt with by numerous Monte Carlo 
integration methods based on efficient sampling of trajectories. The second is 
usually addressed by the use of better discretizations, i.e by substituting better 
effective actions for the naively discetized action. The underlying idea here is 
to construct and use effective actions that lead to improved convergence of 
physical expressions to the (same) continuum limit. 

Typically, physical expression calculated using the naively discretized action 
converge to the continuum as A review of a variety of effective actions 
constructed by improving short-time propagation and using generalizations of 
the Trotter-Suzuki formula [T3|, together with their ranges of application and 
dependence on ordering prescriptions, are given in [T3]. For a long time now, 
the state of the art result has been the convergence of discretized parti- 
tion functions obtained by Takahashi and Imada [15], and Li and Broughton 
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|16] . In these papers the authors used a generahzed form [T7] of the Trotter for- 
mula and the cychc property of the trace to increase the speed of convergence 
of the discretized partition functions to We stress that this increase in 

the speed of convergence only holds for partition functions and not for am- 
plitudes. Some recent results regarding the efficient implementation of PIMC 
algorithms may be found in refs. p^l9ll20|21|22j . In general, the main feature 
of the improved effective action - faster approach to the continuum limit, i.e 
with smaller number of time slices - translates into numerical speedup, simpler 
integration and smaller statistical fluctuations of the quantities calculated. 

In a recent paper [23], we have introduced the concept of the ideal effective 
action S^f. Discretized amplitudes calculated with the ideal effective action do 
not depend on the discretization coarseness N, i.e. they are all equal to the 
sought-after continuum amplitude 

A%{a,b;T) = A{a,b;T) . (2) 

Note that for free particles is in fact just the naively discretized action. 
In [21] we derived the integral equation for the ideal action for a single par- 
ticle moving in one spatial dimension in a general potential. We then solved 
this integral equation as an asymptotic expansion in e. The obtained trunca- 
tion of the ideal effective action S*^^ to order O(e^) systematically improves 
convergence to the continuum to 

A%\a, h- T) = A{a, b; T) + 0(l/iV^') , (3) 

where p = 1 corresponds to the naive action. Up to now the work has focused 
on one particle theories in one dimension. Within this set of theories the out- 
hned procedure has been used to determine explicit expressions for the set of 
effective discretized actions up to level p = 12. The ensuing speedup of several 
orders of magnitude has been verified by extensive numerical simulations [25j. 
The speedup holds for the calculation of all path integrals - for transition am- 
plitudes, partition functions, energy expectation values (in combination with 
appropriate energy estimators), as well as for calculations of energy levels. 
In this paper we extend the above method to the physically relevant case of 
general many-particle non-relativistic systems in arbitrary number of spatial 
dimensions. This is done within a new analytical approach to the construction 
of the effective actions hierarchy. The new approach gives the same effective 
actions as the the old one when applied to one particle one dimensional theory, 
however its computational complexity is lower, making it possible to determine 
effective actions at higher p levels. In particular, the new approach allows for a 
relatively straight-forward extension of the formalism to many-body theories 
and arbitrary dimensions. 
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The present paper gives the analytical derivation of these results, as well as 
a series of Monte Carlo simulations implementing the newly derived effective 
actions and explicitly displaying that the derived speedup indeed holds. The 
paper is organized as follows: Section [2] introduces the hierarchy of effective 
discretized actions and gives the new and extended (many particles, higher 
dimensions) analytical approach to the construction of these effective actions. 
Section [3] presents numerical results that demonstrate the validity of the ana- 
lytically derived speedup in convergence. The efficiency of the new approach is 
further demonstrated in Section H] through the calculation of energy spectra. 
A construction of virial energy estimators corresponding to newly derived ef- 
fective actions and the results of numerical simulations that implement them 
are given in Section [51 

We conclude the paper with a brief summary of obtained results, indicat- 
ing what we see to be the next steps in our line of research and its future 
applications. Appendix A gives a list of integration formulas needed for cal- 
culations up to level p = 5. Appendix B lists the corresponding effective 
actions representing the new state-of-the-art for path integral calculations of 
non-relativistic many-body systems in arbitrary number of dimensions. Higher 
level effective actions can be found on our web site 1261. 



2 New approach to the derivation of effective actions 

We present a method for systematically increasing the efficiency of PIMC cal- 
culations of a non-relativistic quantum system consisting of M distinguishable 
particles in d spatial dimensions, with a Hamiltonian of the form H = K + V , 
where K is the usual kinetic energy operator, and V is the potential describing 
inter-particle interactions and interactions with external fields. The above set 
of theories encompasses a large number of physically relevant models. However, 
as we shall see from the derivations, the method is in principle applicable to all 
quantum theories. For the considered class of theories, the naively discretized 
action Sn is 



Vectors with one index (e.g. g„) represent the set of positions of all particles 
after n time steps of length e = T/N, while vectors with two indices (e.g. 
Qn/) represent dimensional positions of particle i at time step ne. We have 
also introduced the associated discretized velocities 6n/ = Qn+i/ — Ini and 
mid-point coordinates g„ = (g^ -|- g„+i)/2. To summarize, n counts the time 
steps, £ the different particles. In the most compact notation, we may think of 




(4) 
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a configuration (trajectory) of the quantum system as a single vector q wliose 
individual components take on Md possible values, representing positions 
of all the particles at a given time step. In this notation the N ^ oo limit of 
the above discretized amplitude is symbolically written as the path integral 

g(r)=6 

A{a, h-T)= J [dq] e-^[^Wl . (5) 

q{0)=a 



Note that we are using the mid-point ordering prescription and units in which 
h and the particle masses have been set to unity. 

The defining relation for path integrals as the continuum limit of discretized 
amplitudes given by equation ([T]) follows from the completeness relation (de- 
composition of unity) 

A{a, b;T) = / dqi--- dqN-i A{a, qi;e)--- A{qN^i, b; e) , (6) 



through the substitution of short-time amplitudes A(g„, e) calculated to 
first order in time step e, leading to the naive action in equation ([1]). This 
is what gives the convergence to standard path integral expressions. A 
faster converging result may be obtained by evaluating the amplitudes under 
the integral to higher orders in e. From the above completeness relation, it 
follows that the ideal discretized action leads to exact propagation in time, 
and is given in terms of the exact amplitude, according to 

A{qn,qn+i;e) = {2ne) ^ e"*" . (7) 
The ideal discretized action is simply the sum of expressions S*: 

S*N = ^ Sn{(ln,qn+i;£) ■ (8) 
n=0 



Full knowledge of the ideal action S^j is equivalent to knowing the exact ex- 
pression for all the amplitudes A. At first, this would seem to indicate that 
nothing new is to be gained by equation ([7]). This, however, is not the case. We 
will use equation ([7j) to input new analytical information into our numerical 
procedure by calculating amplitudes within some available analytical approx- 
imation scheme. We focus on the calculation of the short time propagation as 
a power series expansion in e, which starts from the naive action (jl]). The de- 
tails of this calculation have been inspired by a similar derivation given in 
One can also attack the problem of solving equation ([7]) using various other 
approximative schemes, especially in the case when short time approximation 
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is not appropriate. The application of the Feynman and Kleinert variational 
approach [27f28| . further developed in [29]|30f31f32ll33] . could prove useful is 
such cases. 

Improved 1 / A^^ convergence of path integrals follows from calculating the gen- 
eral short-time amplitude A{qn, qn+i] £) up to the order e^. To this end, we 
first perform a simple shift of integration variable q = + x about a fixed 
referent trajectory ^, 

x{e/2)=0 

A{q^, e) = e'^-Kl / [dx] ^-^rM¥'^^^-^0) (g) 

a:(-e/2)=0 

We have also shifted the time from t G [ne, {n + l)^] to s G [— e/2,5/2]. The 
referent trajectory ^ satisfies the same boundary conditions as q. As a result, 
the new integration variable x vanishes at the boundaries. The action Sn[i] is 
defined as 

e/2 

Snm= J ds(^^e+vio) , (10) 

-e/2 

and U{x; ^) = V{C, + x) — V{^) — x^, with dots representing derivatives over 
time s. The amplitude may now be written as 

A{qn,qn+i;e) = —-jTj{e^~^^^ ' ''M , (11) 
{2iTe) 2 \ / 

where (...) denotes the expectation value with respect to the free particle 
action. The above expression holds for any choice of referent trajectory C,- 
Equations and ffTTl) now give the expression for the ideal effective action 

S: = S4^]- In (e-G^''''^''-''^y (12) 

The class of theories considered is free of ordering ambiguities, i.e. different 
ordering prescriptions yield the same continuum result. For concreteness we 
work in the mid-point prescription, calculating the above amplitude as a power 
series in time step e. The free particle expectation value in equations ffTTl) and 
( fT2i) is calculated using a Taylor expansion in powers of U: 
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+ -I I dsds' ([/(x;Of/(x';n) + 



(13) 



with shorthand notation x' = x{s'),C,' = C,{s'). By expanding U{x;^) around 
the referent trajectory ^, we get 

U{x- = X, {d.ViO - 6) + lx,x,d,d,V{0 + ■■■ (14) 

From now on we assume summation over repeated indices. The expectation 
values of products (xj(s) . . .Xj{s')) are now calculated in the standard way 
with the help of the free-particle generating functional given in terms of the 
propagator: 

A(.,.')., = ^^(^-^') (|-^) {l + s') 

+ f ^(^'-^)(| + ^)(|-^')- (15) 

From this follow the usual Wick's theorem results {xi{s)) = 0, {xi{s)xj{s')) = 
A{s,s')ij, etc. Note that the generating functional (and so all the expecta- 
tion values) is independent of the specific choice of referent trajectory, i.e. the 
boundary conditions for the x are the same for all choices of ^, and so the 
propagator is always given by ([T^ . Different choices of ^ simplify different 
approximation schemes: using the classical trajectory for ^ is optimal for re- 
covering semi-classical expansion, the choice of linear referent trajectories will 
turn out to be optimal for short-time expansion. 

We next wish to perform the remaining integrations over s. Because of the ex- 
plicit dependence of the referent trajectory on s, we first expand the potential 
and all its derivatives in (fT^ around some reference point. The choice of qn 
for that point corresponds to the mid-point ordering prescription. Once one 
chooses the referent trajectory ^(s), all expectation values in (fTTj) are given in 
terms of quadratures. The choice of linear referent trajectories ^(s) = qn + 
makes all of these integrals solvable in closed form, allowing us to determine 
the ideal effective action. 

Before doing these explicit calculations, we first look at which terms need to 
be retained in order to get the sought-after l/N^ convergence. It is easy to 
show that the ideal effective action is a sum of terms of the form 

e-5^''[d^'Wiqn))---{d^''Viqn)), (16) 
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where a, {3, and 7i,...,7fc arc nonnegative whole numbers constrained by 
simple dimensional analysis to satisfy 

1 ^ 

a + (3^k + -Y^^r. (17) 

^ r=l 



Short time propagation satisfies the diffusion relation S'^ oc s. Thus, l/N^ 
convergence follows from keeping all the terms satisfying a + P < p + 1. An 
equivalent, but practically more useful, form of the criterion determining the 
relevant terms at level p is 

1 ^ 

^+oE7r<P+l. (18) 

^ r=l 



We now continue with the explicit evaluation of the ideal effective action, 
illustrating the general procedure on calculations to level p — 2. The action is 
now 

Sn[^] + n^n) + ^aJy(g-„)) + 0{e') . (19) 



Note that is the i-th component of vector 5„ while dfj is shorthand for 
didj. The first two terms in the above expression correspond to the naive 
action, while the third term gives contributions of order e^. The remaining 
contribution at level p — 2 comes from the expectation value 



2 

e 

2 

_£ 

2 

I ds - {xiXj) didjV{0 + O{e')^ 



= 1 



2 

2 



2 

2 



l-^9V(g„) I dsA{s,s)i,i +0{e'). 



As promised, the remaining integral is easily evaluated in closed form to yield 
£^/ 6. Appendix A lists all the related integrals needed for higher level calcu- 
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lations. The expectation value now equals 



^-Jds uix;0^ ^ ^ _ £_5V(g„) + 0(£3) = e-i^^'^(''") + Oie'') . (20) 



Finally, the p = 2 level discretized effective action is simply 



-(P=2) 
'at 



n=0 



■ (21) 



One similarly obtains the higher p level effective actions. Appendix B gives the 
explicit analytical expressions for many-particle discretized effective actions 
at level p < 5. There are no obstacles in going to higher values of p. The 
derived expressions become algebraically more complex, and calculations are 
best done using a package for symbolic calculus such as Mathematica. Higher 
level effective actions can be found on our web site 1261. 



3 Numerical results 



In this section we present results of numerical PIMC simulations that confirm 
the analytically derived speedup in convergence of discretized path integrals. 
To do this, we have conducted a series of PIMC simulations of transition 
amplitudes for a two-dimensional system of two particles interacting through 
potential 

Vin, r,) = l(ri - + ^in - r,f + |(fi + r,f . (22) 



Numerical simulations, based on our SPEEDUP [26j PIMC code, have been 
performed for different values of couplings gi and g2 and for variety of initial 
and final states. The associated continuum limit amplitudes A^'^'^ have been es- 
timated by fitting polynomials in to the discretized values A^^\ according 
to the analytically derived relation ([3]): 

= A^P^ + ^ + ^ + .... (23) 



For all values of p the fitted continuum values A^^^ agree within the error bars. 
The obtained dependence gives explicit verification of the analytically 
derived increase in convergence. The typical case is illustrated in Figure [H 
The expected convergence may be most easily discerned from Figure 
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Fig. 1. Convergence of discretized amplitudes A^^^ to the continuum as functions of 
N for p = 1,2, 3, 4, 5 for a system of two particles in two dimensions moving in the 
quartic potential given in (|22p with coupling gi = 10, g2 = 0, time of propagation 
T = 1, and initial and final states a = (0, 0; 0.2, 0.5), b = (1, 1; 0.3, 0.6). The number 
of MC samples was 10^. The horizontal dashed line represents the continuum limit, 
sohd lines correspond to the fitted functions (j23p . 

0.01 p ~ ~ ~ ~ ~ ~ ~ ~ n 





Fig. 2. Deviations from the continuum limit — ^| as functions of for p = 1, 3, 5 
for the system of two particles in two dimensions in the quartic potential given in 
()22p with gi = 10, g2 = 0, time of propagation T = 1, initial and final states 
a = (0,0; 0.2, 0.5), b = (1, 1; 0.3, 0.6). The number of MC samples was lO^^ {p = 1), 
10^ (p = 3), 10^^ {p = 5). Solid lines give the leading 1/N'P behavior, dashed lines 
correspond to the fitted functions (p3]) . 
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a plot of the deviations of discretized amplitudes from the continuum limit. 
As can be seen, the increase of level p leads to an ever faster approach to the 
continuum. 

As a result of the newly presented method, the usual simulations (in which 
one calculates specific physical quantities such as the one in Figured]) proceed 
much faster than by using standard methods. On the other hand. Figure [2] 
is itself time consuming since it illustrates subdominant behavior. For this 
reason the figures contain only results obtained by effective actions to level 
p = 5. Note, however, that the p = 5 curve corresponds to a precision of four 
decimal places even for an extremely coarse discretization such as = 2. 

In usual PIMC calculations one always chooses the number of MC samples so 
that the stochastic error of numerical results is of the order of deviations from 
the continuum limit. In Figure [2] the number of MC samples had to be much 
larger, in order for deviations from the continuum limit to be clearly visible. 



4 Energy spectra 

The improved convergence of path integral expressions for amplitudes leads 
directly to the same kind of improvement in the convergence of discretized 
partition functions, owing to the relation 



where the inverse temperature /3 plays the role of the time of propagation T. 
From the previous relation we directly obtain the path-integral presentation 
of the partition function: 



The partition function is the central object for obtaining information about 
statistical systems. In addition, the partition function offers a straightforward 
way for extracting information about low lying energy levels of a system. This 
follows from evaluating the trace in the definition of the partition function in 
the energy eigen-basis 




(24) 




(25) 



oo 



Z{(3) = J2 due 



-f3E, 



(26) 



n=0 
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where En and dn denote corresponding energy levels and degeneracies. A de- 
tailed description of the procedure for extracting energy levels from the par- 
tition function is given in The free energy of the system, given by 

F(/3) = -ilnZ(/3), (27) 



tends to the ground state energy Eq in the large (3 limit. Similarly, we introduce 
auxiliary functions 

i,„£M^sSA£:!i. ,28) 

P dn 



which can be fitted for large /3 to 

=K-^ln(l + ae-'3^). (29) 



and which tend to the corresponding energy level E^- 

In this way, by studying the large f3 behavior of the free energy and the 
corresponding auxiliary functions, one obtains the low lying energy spectrum 
of the model. However, in numerical simulations we are inevitably limited 
to a finite range of inverse temperatures. In addition, the above procedure 
for the construction of the auxiliary functions /„ is recursive, i.e. in order to 
construct we need to know all the energy levels below En- This leads to the 
accumulation of errors as n increases, practically limiting the number of energy 
levels that can be calculated. Note that the orders of magnitude increase in 
precision of the presented method reduces the magnitude of the accumulated 
error and thus allows us to extract a larger number of energy levels. 

Figures [3] and H] illustrate the calculation of low lying energy levels using 
the more efficient PIMC formalism presented in this paper on the case of a 
two-dimensional system of two distinguishable particles interacting through a 
quartic potential fl22|) . Figure [3] demonstrates that the improved convergence 
of free energies is the same as in the case of amplitudes. Table [1] gives the 
calculated energy levels for quartic coupling from gi = (free theory) to 
Qi = 10 (strongly interacting theory). 

Figure H] presents a comparison of ground energy values calculated using the 
derived effective actions with those calculated by perturbative expansion. The 
graph gives the dependence of obtained values of ground energy on the cou- 
pling constant gi. The agreement with perturbative results is excellent for 
small values of the coupling. Around gi ^ 1 even the higher order perturba- 
tive results start deviating substantially from the exact result, while for larger 
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Fig. 3. Convergence of discretized free energies F^' {(3) to the continuum as functions 
of for p = 1, 2, 3, 4, 5 for the system of two particles in two dimensions in potential 
(|22|) for gi = 1, g2 = 1, P = 1. The number of MC samples was 10''. The horizontal 
dashed line represents the exact free energy. 
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Fig. 4. Ground state energy versus coupling constant gi for two-dimensional two- 
particle system in the quartic potential with g2 = 1/9. Numerical simulations were 
performed with 10^ MC samples, level p = 5 effective action and N = 64. Full lines 
give perturbative expansion results of the first order (green), second order (blue), 
and third order (pink). 

couplings we are well in the non-perturbative regime in which such expansions 
are useless. Throughout, the calculations based on the use of the derived ef- 
fective actions converge to the exact ground energy even for relatively rough 



13 



(7i 

yi 


Eq 








n n 


1 .OOO 1 \^ J 


9 ^^^71 ((\\ 

L.OO I 1 \ \j J 






0.1 
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2.497(3) 


2.94(3) 
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2.6327(6) 


3.098(4) 


3.57(3) 





Table 1 

Low lying energy levels for two-dimensional two-particle system interacting through 
potential (j22p with g2 = 1/9. Calculations were done with 10^ MC samples, level 
p = 5 effective action and = 64. The degeneracies of the calculated energy levels 
were found to be do = 1, di = 2, d2 = 3, ds = 6. 

discretizations. From Tabled] we see that this holds even in the case of strongly 
interacting theory. 



5 Energy estimators 

One of the basic properties of physical systems is the internal energy given as 
the energy expectation value. Various ways of evaluating energy expectation 
values in PIMC simulations are based on different types of energy estimators 
- expressions whose thermal expectation values are used to calculate the dis- 
cretized internal energy. The practical choice of the best estimator usually 
depends on the model and on the specific values of physical parameters. It is 
determined through a comparison of the properties of different estimators (e.g. 
variance, convergence, etc.), both in terms of calculated numerical values and 
their functional dependencies. Details of the derivation and characteristics of 
several energy estimators are given in [12] and [35]. The virial energy estimator 
turns out to be most advantageous due to its roughly constant Monte Carlo 
variance with the increase in the number of time slices. In this section we 
construct a hierarchy of virial energy estimators with improved convergence 
to the continuum limit as while keeping the feature of constant variance 

of the naive virial estimator. 

First we briefly review the standard derivation of the virial estimator ^35|36j- 
Starting from the formula U{P) = —dplogZ^P), the straightforward general- 
ization to discretized expressions reads: 

UMP) ^ (30) 

The above partial derivative can be rewritten in terms of e as 9/? = 1/N d^. 
In order to simplify the calculation of this derivative, we remove the e depen- 



14 



dence of the path integral measure and of the kinetic term of the discretized 
expression for by simply rescaling q (l\fs. Now the differentiation over 
e only affects the rescaled potential term in the exponent of the expression for 
Zat. After the differentiation, we reinstall the original variables g, and obtain: 

Vn{^) = (^) ' Jdqi... dqM E^iriai e-^^. (31) 



Here, E^riai is the standard virial energy estimator given as 



1 ^V.,, . 1 



n=0 



Evirial = J^Yl [^iln) + - Qn,! diV{qn) ) • (32) 



This estimator yields a typical 1/A^ convergence. 

As shown in [36], improved convergence of expectation values follows from 
using the appropriate effective actions and energy estimators. As we have 
seen, virial estimators follow directly from the form of the discretized action. 
The level p estimator is obtained by substituting S^^ for Sn in this procedure. 
Writing the ideal effective action and its associated virial energy estimator as 



oo 



S% = S^ + Y.a^^^ (33) 

p=2 

oo 

where cx^^-' and e'-^-' are corresponding contributions proportional to e^. Each 
a^P^ is a sum over all the time-slices, i.e. a^^^ = J2n=o'^n^- The same relation 
holds between e^^^ and e!^\ The outlined procedure now gives the following 
simple connection between ideal action and estimator: 

e(f) = i . (35) 



The analytically derived improvement in convergence has been verified by a 
series of Monte Carlo simulations that have been performed on a two particle 
two dimensional system in a quartic potential (1221) for a variety of different 
values of coupling constant and inverse temperature. Typical results are pre- 
sented in Figures [5] and [61 Figure [5] shows convergence of discretized energy 
expectation values for different levels p. As level p increases, we see that the 
continuum limit is approached faster, with ever smaller values of N. Numerical 
results conform precisely to the analytically derived l/N^ increase in conver- 
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Fig. 5. Convergence of discretized energy expectation values {(5) to the contin- 
uum as functions of for p = 1,2,3,4,5 for the system of two particles in two 
dimensions in a potential (j22p . with gi = 1, g2 = 1/9, /3 = 1. The number of MC 
samples was 10^. The horizontal dashed line represents the exact internal energy. 




Fig. 6. Deviations from the continuum limit (P) — U{P)\ as functions of for 
p = 1,2, 3, 4, 5 for the system of two particles in two dimensions in quartic potential 
(1221) . with gi = I, g2 = 1/9, (3 = 1. The number of MC samples was lO'^ {p = 1), 10^ 
(p = 2), lO^'' {p = 3), 10^^ {p = 4,5). Solid lines give the leading l/N^ behavior, 
dashed lines correspond to the fitted functions (|23p. 
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gence as can be seen from the log-log plot of the deviations of discretized 
values from the continuum limit shown in Figure [H 
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6 Conclusions 

We have presented a derivation of discretized effective actions and energy es- 
timators that lead to substantial, systematic speedup of numerical procedures 
for the calculation of path integrals of a generic many-particle non-relativistic 
theory. The derived speedup holds for all path integrals - for transition am- 
plitudes, partition functions, expectation values, as well as for calculations of 
energy levels. The obtained analytical results have been numerically verified 
through simulations of path integrals for multi-particle model with quartic 
coupling. In the paper we present explicit expressions for the effective actions 
up to level p = 5. The developed calculation scheme has been completed and is 
ready for applications to relevant problems in condensed matter physics. Fur- 
ther analytical work will focus on the generalization of the outlined scheme to 
more complex bosonic and fermionic systems, e.g. field theories. 
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Appendix A: Integrals of propagators 



In this appendix we present the values of all integrals necessary for calculations 
of effective actions at levels p < 5. 



2 

ds A(s,s)^ s"-^'^ 



n odd 

2-l-n ^1-k+n ^(^^ ^ ]-2ttn-^ ^ ^^^^ 



where B is the Euler beta function. All the multiple integrals needed have the 
same boundaries of integration as the preceding case. 

/ dsds' A{s, s') = ^ J dsds' A(s, s'f = ^ 

J dsds' A{s, s')A{s', s') = ^ J dsds' A{s, s')A{s, s') = |^ 

/ dsds' A(s, s') = ^ J dsds' ss' A{s, s') = ^ 

/ dsds' ss''^ A{s, s') = / dsds' s'^s''^ A(s, s') = 

J dsds' ss' A{s, s'y = J dsds' s^ A(s, s') = 

J dsds' s'2 A{s, s'y =2^0 / dsds' s^ A{s', s')A{s, s') - 

J dsds' s2 A(s, s)A{s, s') = -^ J dsds' ss' A{s, s')A{s', s') - 



2520 J ^c.^^ ^ ^y-^ , ^ j^y-^, " J — ^260 

1680 J ^y^, ^ j^y^ , ^ J — 

J dsds' A{s, s)A{s, s')A{s', s') = ig / dsds' A{s, s)A{s, s'f = ^ 

/ dsds' A{s, s')A{s, s f = ^ J dsds'ds" A{s, s")A{s', s") = g 



Appendix B: Effective actions for levels p < 5 



The ideal effective action for a general non-relativistic multi-particle system 
in the mid-point prescription is given in ( l33l) as a sum of terms o"*-^^ containing 
all contributions of order e^. As we have seen, a^^'^ is the sum over all the 
time-slices, i.e. cr^^^ = J2n=o \ where cr^^^ is the contribution of time-slice n. 
In this appendix we give a list of the explicit expressions for cr^^') up to p = 5 
in shorthand notation in which V = V{qn) and St = Sn,i- The expressions for 
higher levels can be found on our web site 
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+ 13440'^'^'^ ^ 1440 '^^^^'^^^■^ + 53760 "^'^-^^"^ ^ + 
322560 

= ^ — - -^d^vd'^d^v - J^did^vdid^v - 

241920 1680 40320 

5 5 5 

^ ^^^^^ - ^^'^^■'^'^ - fS^'''^'^^'^ - 

+ 1935360 53760 ^ 

40320 mjki 32256 mfci ' 

11612160 ^^^■'^'"^"^ ^ 92897280 ^^J*^'™*" 
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